clc
clear;close all
clear functions
%% 导入数据，生成乱序数组
load F:\matalable文件\毕业设计\数据\PCA处理后数据+特征隐藏后数据\data.mat
iris_data=data_finish;         %导入数据
arry=randperm(10830);
%% 划分数据
input_train = iris_data(arry(1:7581),1:4)';
output_train = iris_data(arry(1:7581),5)';
input_test = iris_data(arry(7582:end),1:4)';
output_test = iris_data(arry(7582:end),5)';

% 节点个数
input_num = 4;                    %4个特征
hidden_num = 7;                   %隐藏层神经元
output_num = 4;                   %输出类别
%% 把输出1维转换为4维[1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1]
output=zeros(7581,4);
for i=1:7581
    switch output_train(i)
        case 1
        output(i,:)=[1 0 0 0];
        case 2
        output(i,:)=[0 1 0 0];
        case 3
        output(i,:)=[0 0 1 0];
        case 4
        output(i,:)=[0 0 0 1];
    end
end
[xx,output]=max(output,[],2);

output=categorical(output);
output2=categorical(output_test);
%可以使用ind2vec函数
%% 数据归一化
method=@mapminmax;                 %最大归一化
%method=@mapstd;                   %标准归一化
[input,inputs] = method(input_train);
%[output,outputs] = method(output_train);
input_test_guiyi = method('apply',input_test,inputs);

%%  数据平铺
% 将数据平铺成1维数据只是一种处理方式
% 也可以平铺成2维数据，以及3维数据，需要修改对应模型结构
% 但是应该始终和输入层数据结构保持一致
input =  double(reshape(input, 4, 1, 1, 7581));
input_test_guiyi  =  double(reshape(input_test_guiyi , 4, 1, 1, 3249));
Input=cell(7581,1);
Input_test_guiyi=cell(3249,1);
for i = 1 : 7581
    Input{i, 1} = input(:, 1, 1, i);
end
for i=1:3249
    Input_test_guiyi{i, 1} = input_test_guiyi( :, 1, 1, i);
end

%% 定义GWO算法参数
% optimize_num=2;             %优化参数个数
SearchAgents_no=40;           % 狼群数量，The number of search agents
Max_iteration=50;             % 最大迭代次数，Maximum numbef of iterations
dim=4;                        % 此例需要优化两个参数，number of your variables
lb=[1e-4 1e-4 10 32];         % 参数取值下界，Upper bound of your parameters%%学习率 L2正则化 最小批次 隐藏层神经元
ub=[0.1 0.1 256 128];         % 参数取值上界, Lower bound of your parameters
Convergence_curve=...
zeros(1,Max_iteration);       % record the optimal of each generation
Convergence_parameter=...     % converge each generation optimal variables of alpha
zeros(Max_iteration,dim);
l=1;                          % Loop counter循环计数器,从1开始
%% 取出精确率和召回率
idex_pr=zeros(1,Max_iteration);                            %每代最优的索引
generation_fitness=zeros(1,SearchAgents_no);               %获取每一代中所有个体的标准差
position_total=zeros(SearchAgents_no*Max_iteration,4);     %取出每代中每个狼的位置，用来确定在不同隐藏神经元个数

%% initialize Alpha,beta,and delta position
Alpha_pos=zeros(1,dim);       % initialize the position of alpha
Alpha_score=inf;              % 初始化Alpha狼的目标函数值，change this to -inf for maximization problems
 
Beta_pos=zeros(1,dim);        % initialize the position of beta 
Beta_score=inf;               % 初始化Beta狼的目标函数值，change this to -inf for maximization problems
 
Delta_pos=zeros(1,dim);       % initialize the position of delta
Delta_score=inf;              % 初始化Delta狼的目标函数值，change this to -inf for maximization problems

%Initialize the positions of search agents
Positions=initialize_position(SearchAgents_no,dim,ub,lb);

%% 模型优化
while l<=Max_iteration                            % 对迭代次数循环
    fprintf('Iteration generation number:%d\n',l)
    for m=1:SearchAgents_no                 %取出迭代过程中出现的所有参数组合
       position_total((l-1)*SearchAgents_no+m,:)=Positions(m,:);%取出所有位置
    end

    for i=1:size(Positions,1)                    % 遍历每个狼
    fprintf('Iteration race number:%d\n',i)
%        % Return back the search agents that go beyond the boundaries of the search space
%        % 若搜索位置超过了搜索空间，需要重新回到搜索空间
%         Flag4ub=Positions(i,:)>ub;
%         Flag4lb=Positions(i,:)<lb;
%        % 若狼的位置在最大值和最小值之间，则位置不需要调整，若超出最大值，则回到最大值边界；
%        % 若超出最小值，则回到最小值边界
%         Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb; % ~表示取反          

       % limit position,if position greter than upper bound that return upper bound
       % limit position,if position less than upper bound that return less bound
       % 直接设置边界法处理边界
       for j=1:size(ub,2)
           if Positions(i,j)>ub(j)
               Positions(i,j)=ub(j);
           elseif Positions(i,j)<lb(j)
               Positions(i,j)=lb(j);
           else 
               Positions(i,j)= Positions(i,j);
           end
       end
       % 随机重置法处理边界
%        for j=1:size(ub,2)
%            if Positions(i,j)>ub(j) || Positions(i,j)<lb(j)
%                Positions(i,j)=rand().*(ub(j)-lb(j))+lb(j);
%            else 
%                Positions(i,j)= Positions(i,j);
%            end
%        end
    
        % 计算适应度函数值,caculate fitness function value
       fitness=RNN_Bpfun(Positions(i,:),Input,output,Input_test_guiyi,output_test,SearchAgents_no,Max_iteration);
       generation_fitness(1,i)=fitness;%获取该代中每个狼的适应度值
       if i==size(Positions,1)
           [~,idex_pr(1,l)]=min(generation_fitness);%每一代中最优狼的索引,fitness最小
       end

        % Update Alpha, Beta, and Delta
        if fitness<Alpha_score                       % 如果目标函数值小于Alpha狼的目标函数值
            Alpha_score=fitness;                     % 则将Alpha狼的目标函数值更新为最优目标函数值，Update alpha
            Alpha_pos=Positions(i,:);                % 同时将Alpha狼的位置更新为最优位置
        end
 
        if fitness>Alpha_score && fitness<Beta_score % 如果目标函数值介于于Alpha狼和Beta狼的目标函数值之间
            Beta_score=fitness;                      % 则将Beta狼的目标函数值更新为最优目标函数值，Update beta
            Beta_pos=Positions(i,:);                 % 同时更新Beta狼的位置
        end
 
        if fitness>Alpha_score && fitness>Beta_score &&...
                fitness<Delta_score                  % 如果目标函数值介于于Beta狼和Delta狼的目标函数值之间
            Delta_score=fitness;                     % 则将Delta狼的目标函数值更新为最优目标函数值，Update delta
            Delta_pos=Positions(i,:);                % 同时更新Delta狼的位置
        end
    end
    Convergence_curve(l)=Alpha_score;                % get each generation optimal fitness of alpha
    Convergence_parameter(l,:)=Alpha_pos;            % get each generation optimal variables of alpha
%     a=2-l*((2)/Max_iteration);                       % 对每一次迭代，计算相应的a值，a decreases linearly fron 2 to 0
    a=2*exp(-0.2*l);                                   %指数衰减y=初始值*exp(-k*t)
 
    % Update the Position of search agents including omegas
    for i=1:size(Positions,1)                        % 遍历每个狼
        for j=1:size(Positions,2)                    % 遍历每个维度
 
            % 包围猎物，位置更新
 
            r1=rand();                               % r1 is a random number in [0,1]
            r2=rand();                               % r2 is a random number in [0,1]
 
            A1=2*a*r1-a;                             % 计算系数A
            C1=2*r2;                                 % 计算系数C
 
            % Alpha狼位置更新
            D_alpha=...
                abs(C1*Alpha_pos(j)-Positions(i,j)); 
            X1=Alpha_pos(j)-A1*D_alpha;             
            
            r1=rand();
            r2=rand();
            
            A2=2*a*r1-a;                             % 计算系数A
            C2=2*r2;                                 % 计算系数C
            
            % Beta狼位置更新
            D_beta=...
                abs(C2*Beta_pos(j)-Positions(i,j));      
            X2=Beta_pos(j)-A2*D_beta;           
            
            r1=rand();
            r2=rand();
 
            A3=2*a*r1-a;                             % 计算系数A
            C3=2*r2;                                 % 计算系数C
 
            % Delta狼位置更新
            D_delta=...
                abs(C3*Delta_pos(j)-Positions(i,j)); 
            X3=Delta_pos(j)-A3*D_delta;                  
 
            % 位置更新
            Positions(i,j)=(X1+X2+X3)/3;            
        end
    end

    %当标准差过低，随机进行变异
    std_fitness=std(generation_fitness);%计算一代中的标准差
    if std_fitness<1e-2                %当标准差小于1e-2时，随机变异一部分位置保持狼群多样性
       variation=randperm(SearchAgents_no,randi(SearchAgents_no/2));%随机生成SearchAgents_no/2个1-40的随机变量
       variation_position=initialize_position(SearchAgents_no,dim,ub,lb);
       for v=1:size(variation,2)%把随机选中的位置清零,赋值变异产生的值
           Positions(variation(v),:)=variation_position(variation(v),:);
       end
    else
       Positions=Positions(:,:);%等于原来的值
    end
    %增加迭代次数
    l=l+1;
end


%% 绘图
figure
plot(0:1:50,[-88.9604 Convergence_curve*100],"r-");
xlabel("迭代次数");
ylabel("fitness");
unique_value=unique(Convergence_curve);      %取出fitness中的唯一值从小到大排列
unique_value=unique_value(end:-1:1);         %从大到小排列
unique_idex=zeros(size(unique_value,2),1);   %设置初始索引
for i=1:size(unique_value,2)                 %取出唯一值出现位置
   unique_idex(i)=find(Convergence_curve==unique_value(i),1,"first");%取出唯一值第一次出现的索引值
end
accurate_idex=zeros(size(unique_value,2),1);   %设置初始索引
for i=1:size(unique_value,2)%取出值第一次出现时的索引
  accurate_idex(i)=find(macro_f1==(-unique_value(i)),1,"first");
end

% 取出准确率
accurate_final=zeros(1,length(Convergence_curve));%初始化一个与fitness相同长度+1的数据
for i=1:size(accurate_idex)%平铺数据，把准确率的位置按照fitness对应的位置输出来
    if i==1
      accurate_final(1,1:(unique_idex(i+1)-1))=accurate(accurate_idex(i));
    elseif i==size(accurate_idex,1)
      accurate_final(1,unique_idex(i):end)=accurate(accurate_idex(i));
    else
      accurate_final(1,unique_idex(i):(unique_idex(i+1)-1))=...
      accurate(accurate_idex(i));
    end
end
accurate_final=accurate_final*100;

% 取出精确率和召回率
precision_final=zeros(4,length(Convergence_curve));%初始化一个与fitness相同长度+1的数据
recall_final=zeros(4,length(Convergence_curve));%初始化一个与fitness相同长度+1的数据
for i=1:size(accurate_idex)%平铺数据，把准确率的位置按照fitness对应的位置输出来
    if i==1%单独处理下界
      for j=1:4
         precision_final(j,1:(unique_idex(i+1)-1))=precision(j,accurate_idex(i));
         recall_final(j,1:(unique_idex(i+1)-1))=recall(j,accurate_idex(i));
      end
    elseif i==size(accurate_idex,1)%单独处理上界
      for j=1:4
         precision_final(j,unique_idex(i):end)=precision(j,accurate_idex(i));
         recall_final(j,unique_idex(i):end)=recall(j,accurate_idex(i));
      end
    else
      for j=1:4
         precision_final(j,unique_idex(i):(unique_idex(i+1)-1))=...
         precision(j,accurate_idex(i));
         recall_final(j,unique_idex(i):(unique_idex(i+1)-1))=...
         recall(j,accurate_idex(i));
      end
    end
end
